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Abstract 

The computation of triangular decompositions involves two fundamental 
operations: polynomial GCDs modulo regular chains and regularity test mod- 
ulo saturated ideals. We propose new algorithms for these core operations 
based on modular methods and fast polynomial arithmetic. We rely on new 
results connecting polynomial subresultants and GCDs modulo regular chains. 
We report on extensive experimentation, comparing our code to pre-existing 
Maple implementations, as well as more optimized Magma functions. In 
most cases, our new code outperforms the other packages by several orders of 
magnitude. 

Keywords: Fast polynomial arithmetic, regular chain, regular GCD, subre- 
sultants, triangular decomposition, polynomial systems. 



1 Introduction 



A triangular decomposition of a set F C k[ list of polynomial systems 

Ti, . . . ,T e , called regular chains (or regular systems) and representing the zero set 
V(F) of F. Each regular chain Tj may encode several irreducible components of 
V(F) provided that those share some properties (same dimension, same free vari- 
ables, . . . ). 
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Triangular decomposition methods are based on a univariate and recursive vision of 
multivariate polynomials. Most of their routines manipulate polynomial remainder 
sequences (PRS). Moreover, these methods are usually "factorization free", which 
explains why two different irreducible components may be represented by the same 
regular chain. An essential routine is then to check whether a hypersurface / = 
contains one of the irreducible components encoded by a regular chain T. This is 
achieved by testing whether the polynomial / is a zero-divisor modulo the so-called 
saturated ideal of T. The univariate vision on regular chains allows to perform this 
regularity test by means of GCD computations. However, since the saturated ideal 
of T may not prime, the concept of a GCD used here is not standard. 

The first formal definition of this type of GCDs was given by Kalkbrener in [14] . 
But in fact, GCDs over non-integral domains were already used in several papers 0, 
\TE\ [T2] since the introduction of the celebrated D5 Principle [7j by Delia Dora, 
Dicrescenzo and Duval. Indeed, this brilliant and simple observation allows one 
to carry out over direct product of fields computations that are usually conducted 
over fields. For instance, computing univariate polynomial GCDs by means of the 
Euclidean Algorithm. 

To define a polynomial GCD of two (or more) polynomials modulo a regular chain 
T, Kalkbrener refers to the irreducible components that T represents. In order to 
improve the practical efficiency of those GCD computations by means of subresul- 
tant techniques, Rioboo and the second author proposed a more abstract definition 
in [23]. Their GCD algorithm is, however, limited to regular chains with zero- 
dimensional saturated ideals. 

While Kalkbrener's definition cover the positive dimensional case, his approach can- 
not support triangular decomposition methods solving polynomial systems incre- 
mentally, that is, by solving one equation after another. This is a serious limitation 
since incremental solving is a powerful way to develop efficient sub-algorithms, by 
means of geometrical consideration. The first incremental triangular decomposition 
method was proposed by Lazard in [T5], without proof nor a GCD definition. An- 
other such method was established by the second author in [22] together with a 
formal notion of GCD adapted to the needs of incremental solving. This concept, 
called regular GCD, is reviewed in Section [2J in the context of regular chains. A 
more abstract definition follows. 

Let B be a commutative ring with unity. Let P,Q,G be non-zero univariate poly- 
nomials in M[y]. We say that G is a regular GCD of P,Q if the following three 
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conditions hold: 

(i) the leading coefficient of G is a regular element of B, 
(ii) G lies in the ideal generated by P and Q in M[y], and 

(in) if G has positive degree w.r.t. y, then G pseudo-divides both of P and Q, that 
is, the pseudo-remainders prem(P, G) and prem(Q, G) are null. 

In the context of regular chains, the ring B is the residue class ring of a polynomial 
ring A := k[xi, . . . ,x n ] (over a field k) by the saturated ideal sat(T) of a regular 
chain T. Even if the leading coefficients of P,Q are regular and sat(T) is radical, 
the polynomials P,Q may not necessarily admit a regular GCD (unless sat(T) is 
prime). However, by splitting T into several regular chains Ti, . . . ,T e (in a sense 
specified in Section [2]) one can compute a regular GCD of P, Q over each of the ring 
A/sat (Ti), as shown in [22] . 

In this paper, we propose a new algorithm for this task, together with a theoretical 
study and implementation report, providing dramatic improvements w.r.t. previous 
work [TJJ[22]. Section [3] exhibits sufficient conditions for a subresultant polynomial 
of P, Q E A[y] (regarded as univariate polynomials in y) to be a regular GCD of 
P, Q w.r.t. T. Some of these properties could be known, but we could not find a 
reference for them, in particular when sat(T) is not radical. These results reduce 
the computation of regular GCDs to that of subresultant chains, see Section @] for 
details. 

Since Euclidean-like algorithms tend to densify computations, we consider an eval- 
uation/interpolation scheme based on FFT techniques for computing subresultant 
chains. In addition, we observe that, while computing triangular decomposition, 
whenever a regular GCD of P and Q w.r.t. T is needed, the resultant of P and 
Q w.r.t. y is likely to be computed too. This suggests to organize calculations in 
a way that the subresultant chain of P and Q is computed only once. Moreover, 
we wish to follow a successful principle introduced in [20]: compute in k[xi, . . . , x n ] 
instead of k[xi, . . . , x n ]/sat(T), as much as possible, while controlling expression 
swell. These three requirements targeting efficiency are satisfied by the implemen- 
tation techniques of Section 15.11 The use of fast arithmetic for computing regular 
GCDs was proposed in [6] for regular chains with zero- dimensional radical saturated 
ideals. However this method does not meet our other two requirements and does not 
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apply to arbitrary regular chains. We state complexity results for the algorithms of 
this paper in Sections 15.11 and 15.21 

Efficient implementation is the main objective of our work. We explain in Section I5~51 
how we create opportunities for using modular methods and fast arithmetic in oper- 
ations modulo regular chains, such as regular GCD computation and regularity test. 
The experimental results of Section [6] illustrate the high efficiency of our algorithms. 
We obtain speed-up factors of several orders of magnitude w.r.t. the algorithms 
of [22] for regular GCD computations and regularity test. Our code compares and 
often outperforms packages with similar specifications in Maple and Magma. 



2 Preliminaries 

Let k be a field and let k[x] = k[xi, . . . , x n ] be the ring of polynomials with coeffi- 
cients in k, with ordered variables xx -< ■ • ■ -< x n . Let k be the algebraic closure of 
k. If u is a subset of x then k(u) denotes the fraction field of k[u]. For F C k[x], 
we denote by (F) the ideal it generates in k[x] and by yjF) the radical of (F). For 
H G k[x], the saturated ideal of (F) w.r.t. H, denoted by (F) : H°°, is the ideal 
{Q G k[x] | 3m G N s.t. H m Q G (F)}. A polynomial P G k[x] is a zero-divisor 
modulo (F) if there exists a polynomial Q such that PQ G (F), and neither P nor 
Q belongs to (F). The polynomial P is regular modulo (F) if it is neither zero, nor 
a zero-divisor modulo (F). We denote by V(F) the zero set (or algebraic variety) of 
F in k . For a subset Wck , we denote by W its closure in the Zariski topology. 

2.1 Regular chains and related notions 

Main variable and initial. If P G k[x] is a non-constant polynomial, the largest 
variable appearing in P is called the main variable of P and is denoted by mvar(P). 
The leading coefficient of P w.r.t. mvar(P) is its initial, written init(P) whereas 
lc(P, v) is the leading coefficient of P w.r.t. v G x. 

Triangular Set. A subset T of non-constant polynomials of k[x] is a triangular set if 
the polynomials in T have pairwise distinct main variables. Denote by mvar(T) the 
set of all mvar(P) for P G T. A variable v G x is algebraic w.r.t. T if v G mvar(T); 
otherwise it is free. For a variable v G x we denote by T <v (resp. T >v ) the subsets 
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of T consisting of the polynomials with main variable less than (resp. greater than) 
v. If v G mvar(T), we denote by T v the polynomial P G T with main variable v. 
For T not empty, T max denotes the polynomial of T with largest main variable. 

Quasi-component and saturated ideal. Given a triangular set T in k[x], denote 
by fix the product of the init(P) for all PeT. The quasi- component W(T) of T is 
V(T) \ V(Ht), that is, the set of the points of V(T) which do not cancel any of the 
initials of T. We denote by sat(T) the saturated ideal of T, defined as follows: if T 
is empty then sat(T) is the trivial ideal (0); otherwise it is the ideal (T) : h™- 

Regular chain. A triangular set T is a regular chain if either T is empty, or 
T \ {T max } is a regular chain and the initial of T max is regular with respect to 
sat(T \ {T max }). In this latter case, sat(T) is a proper ideal of k[x]. From now on 
T C k[x] is a regular chain; moreover we write m — |T|, 6 = mvar(T) and u = x\fi. 
The ideal sat(T) enjoys several properties. First, its zero-set equals W(T). Second, 
the ideal sat(T) is unmixed with dimension n — m. Moreover, any prime ideal p asso- 
ciated to sat(T) satisfies p nk[u] = (0). Third, if n = m, then sat(T) is simply (T). 
Given P G k[x] the pseudo-remainder (resp. iterated resultant) of P w.r.t. T, de- 
noted by prem(P, T) (resp. res(P, T)) is defined as follows. If P G k or no variables 
of P is algebraic w.r.t. T, then prem(P, T) = P (resp. res(P, T) = P). Otherwise, 
we set prem(P, T) = prem(P, T <v ) (resp. res(P, T) = res(P, T <v )) where v is the 
largest variable of P which is algebraic w.r.t. T and R is the pseudo-remainder 
(resp. resultant) of P and T v w.r.t. v. We have: P is null (resp. regular) w.r.t. 
sat(T) if and only if prem(P, T) = (resp. res(P, T) ^ 0). 

Regular GCD. Let / be the ideal generated by a/ sat (T) in k(u)[B]. Then L(T) := 
k(u)[6]/J is a direct product of fields. It follows that every pair of univariate polyno- 
mials P, Q G L(T)[y] possesses a GCD in the sense of [23]. The following GCD no- 
tion |22j is convenient since it avoids considering radical ideals. Let T C k[ 
be a regular chain and let P, Q G k[x, y] be non-constant polynomials both with main 
variable y. Assume that the initials of P and Q are regular modulo sat(T). A non- 
zero polynomial G G k[x, y] is a regular GCD of P,Q w.r.t. T if these conditions 
hold: 

(i) \c(G,y) is regular with respect to sat(T); 
(ii) there exist u, v G k[x, y] such that g — up — vt G sat(T); 
(Hi) if deg(G, y) > holds, then (P, Q) C sat(T U G). 
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In this case, the polynomial G has several properties. First, it is regular with 
respect to sat(T). Moreover, if sat(T) is radical and deg(G, y) > holds, then the 
ideals (P,Q) and (G) of L(T)[y] are equal, so that G is a GCD of (P,Q) w.r.t. 
T in the sense of [23J. The notion of a regular GCD can be used to compute 
intersections of algebraic varieties. As an example we will use Formula (pQ) which 
follows from Theorem 32 in [22] • Assume that the regular chain T is simply {R} 
where R = res(P, Q, y), for R ^ k, and let H be the product of the initials of P and 
Q. Then, we have: 

V(P, Q) = W(R,G) U V(H, P, Q). (1) 

Splitting. Two polynomials P, Q may not necessarily admit a regular GCD w.r.t. 
a regular chain T, unless sat(T) is prime, see Example 1 in Section [31 However, if T 
"splits" into several regular chains, then P, Q may admit a regular GCD w.r.t. each 
of them. This requires a notation. For non-empty regular chains T, T 1; . . . , T e C 
k[x] we write T — > (Ti, . . . , T e ) whenever ^/sat(T) = ^/sat(Ti) fl • • • fl ^/sat(T e ), 
mvar(T) = mvar(Tj) and sat(T) C sat(Tj) hold for all 1 < % < e. If this holds, 
observe that any polynomial if regular w.r.t sat(T) is also regular w.r.t. sat(Tj) for 
all 1< % < e. 



2.2 Fundamental operations on regular chains 

We list below the specifications of the fundamental operations on regular chains 
used in this paper. The names and specifications of these operations are the same 
as in the RegularChains library [T5] in MAPLE. 

Regularize. For a regular chain T C k[x] and P in k[x], the operation Regularize(P, T) 
returns regular chains Ti, . . . ,T e of k[x] such that, for each 1 < i < e, P is either 
zero or regular modulo sat(Tj) and we have T — ^(Ti, ■ • • , T e ). 

RegularGcd. Let T be a regular chain and let P,Q G k[x, y] be non-constant 
with mvar(P) = mvar(Q) ^ mvar(T) and such that both init(P) and init(Q) are 
regular w.r.t. sat(T). Then, the operation RegularGcd (P, Q, T) returns a sequence 
(Gi, Ti), . . . , (G e , T e ), called a regular GCD sequence, where Gi, . . . , G e are polyno- 
mials and Ti, . . . , T e are regular chains of k[x], such that T — ►(Ti, . . . , T e ) holds and 
Gi is a regular GCD of P, Q w.r.t. Tj for all 1 < % < e. 
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NormalForm. Let T be a zero-dimensional normalized regular chain, that is, a 
regular chain whose saturated ideal is zero-dimensional and whose initials are all 
in the base field k. Observe that T is a lexicographic Grobner basis. Then, for 
P G k[x], the operation Norma I Form (P, T) returns the normal form of P w.r.t. T in 
the sense of Grobner bases. 

Normalize. Let T be a regular chain such that each variable occurring in T belongs 
to mvar(T). Let P 6 k[x] be non-constant with initial H regular w.r.t. (T). 
Assume each variable of H belongs to mvar(T). Then H is invertible modulo (T) 
and Normalize(P, T) returns NormalForm(iJ _1 P, T) where H~ l is the inverse of H 
modulo (T). 

2.3 Subresultants 

We follow the presentation of [8], [25] and [ID] . 

Determinantal polynomial. Let B be a commutative ring with identity and let 
m < n be positive integers. Let M be a m x n matrix with coefficients in B. Let Mj 
be the square submatrix of M consisting of the first m — 1 columns of M and the 
i-th column of M, for i — m • ■ • n; let det Mj be the determinant of Mj. We denote 
by dpol(M) the element of B[y], called the determinantal polynomial of M, given by 

det M m y n - m + det M m+1 y n_m_1 + • ■ • + det M n . 

Note that if dpol(M) is not zero then its degree is at most n — m. Let Pi, ... , P m be 
polynomials of M[y] of degree less than n. We denote by mat (Pi, . . . , P m ) the m x n 
matrix whose z-th row contains the coefficients of Pj, sorting in order of decreasing 
degree, and such that P« is treated as a polynomial of degree n — 1. We denote by 
dpol(Pi, . . . , P m ) the determinantal polynomial of mat (Pi, . . . , P m ). 

Subresultant. Let P,Q G B[?/] be non-constant polynomials of respective degrees 
p, q with q < p. Let d be an integer with < d < q. Then the d-ih subresultant of 
P and Q, denoted by Sd(P, Q), is 

dpol ( yg -d-l P) y1-d-2p^ . . . ; P) yP-d-lQ, ...,Q). 

This is a polynomial which belongs to the ideal generated by P and Q in B[y]. In 
particular, S (P, Q) is res(P, Q), the resultant of P and Q. Observe that if Sd(P, Q) 
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is not zero then its degree is at most d. When Sd(P, Q) has degree d, it is said 
non-defective or regular, when Sd(P, Q) 7^ and deg (Srf(P, Q)) < d, Sd(P, Q) is said 
defective. We denote by Sd the coefficient of Sd(P, Q) in y d . For convenience, we 
extend the definition to the q-th subresultant as follows: 

g (p q) - I l(Q)Q> if p> q or \c(Q) G B is regular 
9 1 undefined, otherwise 

where 7(Q) = \c(Q) p ~ q ~ . Note that when p equals q and \c(Q) is a regular element 
in B, S q (P, Q) = lc(Q)~ 1 Q is in fact a polynomial over the total fraction ring of B. 

We call specialization property of subresultants the following statement. Let D be 
another commutative ring with identity and \I/ a ring homomorphism from B to D 
such that we have \I/(lc(P)) ^ and ^(lc(Q)) 7^ 0. Then we also have 

S d (*(P),*(Q)) = *(S d (P,Q)). 

Divisibility relations of subresultants. The subresultants S q -i(P, Q), S q -2 (P, Q), 
. . ., Sq(P, Q) satisfy relations which induce an Euclidean-like algorithm for comput- 
ing them. Following [8] we first assume that B is an integral domain. In the above, 
we simply write Sd instead of Sd{P, Q), for d = q — 1, . . . , 0. We write A ~ B for 
A,B G B[y] whenever they are associated over fr(B), the field of fractions of B. For 
d = q — 1,...,1, we have: 

(r q -i) Sq-i = prem(P, —Q), the pseudo-remainder of P by —Q, 
( r <g-i) if <3q-i 7^ 0, with e = deg(S' 9 _i), then the following holds: prem(Q, —S g -i) = 

lc(Q) (P ~ ,)(,-e)+1 5 e _l, 

(r e ) if Sd_i 7^ 0, with e = deg(Sd_i) < d — 1, thus S^-i is defective, and we have 
(i) deg(S'rf) = d, thus S d is non-defective, 

(ii) Sd-i ~ and \c(Sd-i) d e l Sd-i = Sd d ~ e ~ 1 S e , thus S e is non-defective, 
(m) = Sd-3 = ■ ■ ■ = S e+ i = 0, 

(r e _i) if S'rf and S^-i are nonzero, with respective degrees d and e, then we have 
prem(5' rf , -S d -i) = lc(S d ) ~ e+1 S e -i, 
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We consider now the case where B is an arbitrary commutative ring, following 
Theorem 4.3 in [TQ] . If Sd, Sd-i are non zero, with respective degrees d and e and 
if Sd is regular in B then we have lc(Sd-i) ~ e 1 Sd-i = Sd d ~ e ~ 1 S e ; moreover, there 
exists Cd G M[y] such that we have: 

(-lY-HciSd-^Sd + CdSd-i = HSdfS^. 

In addition Sd-2 = S<z-3 = ■ ■ • = S e+ i = also holds. 



3 Regular GCDs 

Throughout this section, we assume n > 1 and we consider P,Qe k[xi, . . . , x n+ i) 
non-constant polynomials with the same main variable y := x n+ \ and such that 
p := deg(P, y) > q := deg(Q, y) holds. We denote by R the resultant of P and Q 
w.r.t. y. Let T C k[xi, . . . ,x n ] be a non-empty regular chain such that R G sat(T) 
and the initials of P, Q are regular w.r.t. sat(T). We denote by A and B the rings 
k[xi,...,x n ] and k[x%, . . . , x n ]/sat(T), respectively. Let ^ be both the canonical 
ring homomorphism from A to B and the ring homomorphism it induces from A[y] 
to M[y]. For < j < q, we denote by Sj the j-th subresultant of P, Q in A[y]. 

Let d be an index in the range 1- • -q such that Sj G sat(T) for all < j < d. 
Lemma [3] and Lemma H] exhibit conditions under which Sd is a regular GCD of P 
and Q w.r.t. T. Lemma [T] and Lemma [2] investigate the properties of Sd when 
\c(Sd,y) is regular modulo sat(T) and lc(Sd,y) G sat(T) respectively. 

Lemma 1 If\c(Sd,y) is regular modulo sat(T), then the polynomial Sd is a non- 
defective subresultant of P and Q over A. Consequently, ^(Sd) is a non-defective 
subresultant of ^ (P) and^(Q) inM[y]. 

Proof. When d = q holds, we are done. Assume d < q. Suppose Sd is defective, 
that is, deg(Sd, y) = e < d. According to item (r e ) in the divisibility relations of 
subresultants, there exists a non-defective subresultant Sd+i such that 

\c(Sd,y) d ~ e S d = s d d ~ e 1 S e , 

where Sd+i is the leading coefficient of Sd+i in y. By our assumptions, S e belongs 
to sat(T), thus lc(Sd,y) d ~ e Sd G sat(T) holds. It follows from the fact \c(Sd,y) is 
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regular modulo sat(T) that Sd is also in sat(T). However the fact that \c(Sd,y) = 
init(Sd) is regular modulo sat(T) also implies that Sd is regular modulo sat(T). A 
contradiction. □ 



Lemma 2 Iflc(Sd,y) is contained in sat(T), then all the coefficients of Sd regarded 
as a univariate polynomial in y are nilpotent modulo sat(T). 

Proof. If the leading coefficient lc(Sd, y) is in sat(T), then lc(Sd, y) G p holds for all 
the associated primes p of sat(T). By the Block Structure Theorem of subresultants 
(Theorem 7.9.1 of [21J) over an integral domain \s.[x%, . . . ,x n -i)/p, Sd must belong 
to p. Hence we have Sd € a/ sat(T). Indeed, in a commutative ring, the radical of 
an ideal equals the intersection of all its associated primes. Thus Sd is nilpotent 
modulo sat(T). It follows from Exercise 2 of [T] that all the coefficients of Sd in y 
are also nilpotent modulo sat(T). □ 

Lemma [2] implies that, whenever lc(Sd,y) £ sat(T) holds, the polynomial Sd will 
vanish on all the components of sat(T) after splitting T sufficiently. This is the key 
reason why Lemma [T] can be applied for computing regular GCDs. Indeed, up to 
splitting via the operation Regularize, one can always assume that either \c(Sd,y) 
is regular modulo sat(T) or \c(Sd,y) belongs to sat(T). Hence, from Lemma [2] and 
up to splitting, one can assume that either lc(Sd,y) is regular modulo sat(T) or Sd 
belongs to sat(T). Therefore, if Sd & sat(T), we consider the subresultant Sd as a 
candidate regular GCD of P and Q modulo sat(T). 

Example 1 Iflc(Sd,y) is not regular modulo sat(T) then Sd may be defective. Con- 
sider for instance the polynomials P = x\x\ — x\ and Q = x\x\ — x\ in Q[x\, a?2, x?f\ . 
We have prem(P, -Q) = {x\ - x 6 2 ) and R = {x\ - x 6 2 ) 2 . Let T = {R}. The last 
subresultant of P,Q modulo sat(T) is prem(P, — Q), which has degree w.r.t x 3 , 
although its index is 1. Note that prem(P, —Q) is nilpotent modulo sat(T). 

In what follows, we give sufficient conditions for the subresultant Sd to be a regular 
GCD of P and Q w.r.t. T. When sat(T) is a radical ideal, Lemma H] states that 
the assumptions of Lemma [1] are sufficient. This lemma validates the search for a 
regular GCD of P and Q w.r.t. T in a bottom-up style, from So up to Se for some I. 
Lemma [3] covers the case where sat(T) is not radical and states that Sd is a regular 
GCD of P and Q modulo T, provided that Sd, satisfies the conditions of Lemma [T] 
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and provided that, for all d < k < q, the coefficient Sk of y k in Sk is either null or 
regular modulo sat(T). 

Lemma 3 We reuse the notations and assumptions of Lemma [H Then Sd is a 
regular GCD of P and Q modulo sat(T), if for all d < k < q, the coefficient Sk of 
y k in Sk is either null or regular modulo sat(T). 

Proof. There are three conditions to satisfy for Sd to be a regular gcd of P and Q 
modulo sat(T): 

(1) \c(Sd) is regular modulo sat(T); 

(2) there exists polynomials u and v such that Sd — uP — vQ G sat(T); and 

(3) both P and Q are in X := sat(T U 

We will prove the lemma in three steps. We write \P(r) as f for brevity 0. 

Claim 1: If Q and S q ~i are in sat(T), then Sd is a regular gcd of P and Q modulo 
sat(T). 

Indeed, the properties of Sd imply Conditions (1) and (2) and we only need to show 
that the Condition (3) also holds. If d = q holds, then S q -i G sat(T) and we are 
done. Otherwise, S q -i = prem(P, —Q) is not null modulo sat(T), because S q -i = 
implies that all subresultants of P and Q with index less than q vanish over B. If 
both S g := Q and S q _i = prem(P, —Q) are in X, then P is also in X, since lc(Q) is 
regular modulo sat(T) and hence is regular modulo X. This completes the proof of 
Claim 1. 

In order to prove that Q and S q -i are in sat(T), we define the following set of indices 

J = {j | d < j < g,coeff(^,^) i sat(T)}. 

By assumption, coeff(5j, y J ) is regular modulo sat(T) for each j 6 J . Our arguments 
rely on the Block Structure Theorem (BST) over an arbitrary ring [10J and Ducos' 
subresultant algorithm [8], [22] along with the specialization property of subresultants. 

1 We note that the degree of Sk may be less than the degree of Sk, since its leading coefficient 
could be in sat(T). Hence, lc(iSfc) may differ from lc(5fe). We carefully distinguish them when the 
leading coefficient of a subresultant is not regular in B. 
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Claim 2: If J = 0, then Si el holds for all d < i < q. 

Indeed, the BST over B implies that there exists at most one subresultant Sj such 
that d < j < q and Sj ^ sat(T). Therefore all but S q -x are in sat(T), and thus S q ~i 
is defective of degree d. More precisely, the BST over B implies 

l^Sq-tYSq^x = \c(Sq) e S d mod sat(T) (2) 

for some integer e > 0. According to Relation ([2]), lc(5 , q _i) is regular in B. Hence, 
we have S q _i G X. From the definition of S d , we have prem(S'q, — 5* g _i, y) G sat(T). 
This implies S q G X. This completes the proof of Claim 2. 

Now we consider the case J ^ 0. Write J explicitly as J = {jo, ji, ■ ■ ■ with 
£ = \ J~\ and we assume jo < ji < • ■ ■ < je-i- For convenience, we write jg := q. For 
each integer k satisfying < k < I we denote by Vk the following property: 

Si G X, for all d < i < jk- 
Claim 3: The property Vk holds for all < k < £. 

We proceed by induction on < k < i. The base case is k = 0. We need to show 
Si G X for all d < i < jo- By the definition of jo, Sj is a non-defective subresultant 
of P and Q, and coeff(S'j, ?/*) is in sat(T) for all d < i < jo- By the BST over B, 
there is at most one d < i < jo such that Si sat(T). If no such a subresultant 
exists, then we know that prem(Sj , — Sd) is in sat(T). Consequently, Sj G X holds, 
which implies Sj G X for all d < i < jo- On the other hand, if Si is not in sat(T) 
for some d <iq < jo, then Si is similar to Sd over B. To be more precise, we have 

\c(S io ) e S io = lc(S jo ) e S d mod sat(T) (3) 

for some integer e > 0. With the same reasoning as in the case J = 0, we know 
that lc(S'j ) is regular modulo sat(T) and we deduce that S io G X holds. Also, we 
have prem(Sj , — Si ) G sat(T), by definition of S d - This implies Sj G X from the 
fact that lc(5i ) is regular modulo sat(T) (and thus regular modulo X). Hence, we 
have Si G X for all d < i < jo, as desired. Therefore the property Vk holds for k = 0. 

Now we assume that the property Vk-i holds for some 1 < k < i. We prove that 
Vk also holds. According to the BST over B, we know that there exists at most one 
subresultant between Sj k _ x and Sj k , both of which are non-defective subresultants 
of P and Q. If Si G sat(T) holds for all jk-i < i < jk, then we have 

prem(S jk , -Sj^) = \c(S jk ) e S u mod sat(T) 
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for some d < u < jk-\ and some integer e > 0. Thus, we have prem(Sj k , —Sj k _ 1 ) G X 
by our induction hypothesis, and consequently, Sj k G X holds. On the other hand, 
if all subresultants Si (for jk-i < i < jk) but Si k (for some index ik such that 
jk-i < ik < jk) are in sat(T), then Si k is similar to over B, that is, we have 

\c(S ik ) e S ik = lc(5 j J e 5 jfc _ 1 mod sat(T) (4) 

for some integer e > 0. By Relation (j4]), lc(Si k ) is regular modulo sat(T), and thus 
is regular modulo X. Using Relation (j3J) again, we have Si k G X, since S , J - fc _ 1 is in X. 
Also, we have 

prem(S jk , -S ik ) = \c(S jk ) e S u mod sat(T) 

for some d < u < jk-i and some integer e > 0. By the induction hypothesis, we 
deduce S u G X, which implies Sj k G X together with the fact that lc(Si k ) is regular 
modulo X. This shows that Si G X holds for all d < i < jk- Therefore, property Vk 
holds. 

Finally, we apply Claim 3 with k = £, leading to Si G X for all d < i < jg = q, which 
completes the proof of our lemma. □ 

The consequence of the above corollary is that we ensure that Sd is a regular gcd 
after checking that the leading coefficients of all non-defective subresultants above 
Sd, are either null or regular modulo sat(T). Therefore, one may be able to conclude 
that Sd is a regular GCD simply after checking the coefficients "along the diagonal" 
of the pictorial representation of the subresultants of P and Q, see Figure [H 

Lemma 4 With the assumptions of LemmaU\ assume sat(T) radical. Then, Sd is 
a regular GCD ofP,Q w.r.t. T. 

Proof. As for Lemma [31 it suffices to check that P and Q belong to sat(T U {Sd})- 
Let p be any prime ideal associated with sat(T). Define D = k[xi, . . . ,y]/p and let 
L be the fraction field of the integral domain D. Clearly Sd is the last subresultant 
of P, Q in B>[y] and thus in L[y]. Hence Sd is a GCD of P, Q in L[y]. Thus Sd 
divides P, Q in L[y] and pseudo-divides P, Q in Therefore prem(P, Sd) and 

prem(Q, Sd) belong to p. Finally prem(P, Sd) and prem(Q, Sd) belong to sat(T). 
Indeed, sat(T) being radical, it is the intersection of its associated primes. □ 
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Figure 1: A possible configuration of the subresultant chain of P and Q. In the left, 
P and Q have five nonzero subresultants over k[x], four of which are non-defective 
and one of which is defective. Let T be a regular chain in k[x] such that lc(P) 
and \c(Q) are regular modulo sat(T). Further, we assume that lc(Si) and \c{S^) 
are regular modulo sat(T), however, lc(S 6 ) is in sat(T). The right hand side is a 
possible configuration of the subresultant chain of P and Q. In the proof of Claim 
3, the set J is {jo = 4} and ji = 7, whereas %q = 2 and i\ — 6 are the indices of 
defective subresultants over k[x]/sat(T). In this case, Si is a regular gcd of P and 
Q modulo sat(T). 
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4 A regular GCD algorithm 



Following the notations and assumptions of Section [3] we propose an algorithm for 
computing a regular GCD sequence of P, Q w.r.t. T. as specified in Section 12.21 
Then, we show how to relax the assumption R G sat(T). 

There are three main ideas behind this algorithm. First, the subresultants of P, Q 
in A[y] are assumed to be known. We explain in Section how we compute them 
in our implementation. Secondly, we rely on the Regularize operation specified in 
Section I2~2l Lastly, we inspect the subresultant chain of P, Q in A[y] in a bottom-up 
manner. Therefore, we view Si, 52,... as successive candidates and apply either 
Lemma HI (if sat(T) is known to be radical) or Lemma [3j 

Case where R G sat(T). By virtue of Lemma [1] and Lemma [2] there exist regular 
chains 7\, . . . , T e C k[x] such that T — ► (Ti, . . . , T e ) holds and for each 1 < i < e 
there exists an index 1 < d{ < q such that the leading coefficient k^S,^,?/) of the 
subresultant S^ is regular modulo sat(T) and Sj G sat(Tj) for all < j < di. Such 
regular chains can be computed using the operation Regularize. If each sat(T) is 
radical then it follows from Lemma |4] that (S^T), . . . , (Sd e ,T e ) is a regular GCD 
sequence of P,Q w.r.t. T. In practice, when sat(T) is radical then so are all 
sat (Tj), see [2]. If some sat(T) is not known to be radical, then one can compute 
regular chains T it x, . . . ,T i>ei C k[x] such that T — > (T^i, . . . ,T je J holds and for 
each 1 < £i < there exists an index 1 < d^ < q such that Lemma [3] applies 
and shows that the subresultant Sd e . is regular GCD of P, Q w.r.t. T^. Such 
computation relies again on Regularize. 

Case where R sat(T). We explain how to relax the assumption R G sat(T) 
and thus obtain a general algorithm for the operation RegularGcd. The principle is 
straightforward. Let R = res(P,Q,y). We call Regularize^, T) obtaining regular 
chains Ti, . . . , T e such that T — > (T\, . . . , T e ). For each 1 < i < e we compute 
a regular GCD sequence of P and Q w.r.t. T as follows: If R G sat(T) holds 
then we proceed as described above; otherwise R G" sat(T) holds and the resultant 
R is actually a regular GCD of P and Q w.r.t. T by definition. Observe that 
when R G sat(T) holds the subresultant chain of P and Q in y is used to compute 
their regular GCD w.r.t. Tj. This is one of the motivations for the implementation 
techniques described in Section [5j 
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5 Implementation and Complexity 



In this section we address implementation techniques and complexity issues. We 
follow the notations introduced in Section [3j However we do not assume that R = 
res(P, Q, y) belongs to the saturated ideal of the regular chain T. 

In Section I5TT1 we describe our encoding of the subresultant chain of P, Q in k[x][y]. 
This representation is used in our implementation and complexity results. For sim- 
plicity our analysis is restricted to the case where k is a finite field whose "charac- 
teristic is large enough". The case where k is the field Q of rational numbers could 
be handled in a similar fashion, with the necessary adjustments. 

One motivation for the design of the techniques presented in this paper is the solv- 
ing of systems of two equations, say P = Q = 0. Indeed, this can be seen as a 
fundamental operation in incremental methods for solving systems of polynomial 
equations, such as the one of [22J. We make two simple observations. Formula [1] 
p. |6] shows that solving this system reduces "essentially" to computing R and a 
regular GCD sequence of P, Q modulo {R}, when R is not constant. This is par- 
ticularly true when n = 2 since in this case the variety V(H, P, Q) is likely to be 
empty for "generic" polynomials P, Q. The second observation is that, under the 
same genericity assumptions, a regular GCD G of P, Q w.r.t. {R} is likely to exist 
and have degree one w.r.t. y. Therefore, once the subresultant chain of P, Q w.r.t. 
y is calculated, one can obtain G "essentially" at no additional cost. Section 15.21 
extends these observations with complexity results. 

In Section an algorithm for Regularize and its implementation are discussed. We 
show how to create opportunities for using fast polynomial arithmetic and modular 
techniques, thus bringing significant improvements w.r.t. other algorithms for the 
same operation, as shown in Section [6j 

5.1 Subresultant chain encoding 

Following we evaluate (xi, . . . , x n ) at sufficiently may points such that the sub- 
resultants of P and Q (regarded as univariate polynomials in y = x n+ \) can be 
computed by interpolation. To be more precise, we need some notations. Let di be 
the maximum of the degrees of P and Q in Xj, for alH = 1, . . . , n + 1. Observe that 
hi := 2did n is an upper bound for the degree of R (or any subresultant of P and Q) 



16 



in Xi, for all i. Let B be the product (&i + 1) • • • (b n + 1). 

We proceed by evaluation / interpolation; our sample points are chosen on an n- 
dimensional rectangular grid. We call "Scube" the values of the subresultant chain 
of P, Q on this grid, which is precisely how the subresultants of P, Q are encoded 
in our implementation. Of course, the validity of this approach requires that our 
evaluation points cancel no initials of P or Q. Even though finding such points 
deterministically is a difficult problem, this created no issue in our implementation. 
Whenever possible (typically, over suitable finite fields), we choose roots of unity 
as sample points, so that we can use FFT (or van der Hoeven's Truncated Fourier 
Transform [13J); otherwise, the standard fast evaluation / interpolation algorithms 
are used. We have 0(d n+1 ) evaluations and 0(d 2 l+1 ) interpolations to perform. Since 
our sample points lie on a grid, the total cost becomes 

oL^pom) or ofB^t^^), 

depending on the choice of the sample points (see e.g. [24) for similar estimates). 
Here, as usual, M(6) stands for the cost of multiplying polynomials of degree less 
than 6, see [UJ Chap. 8]. Using the estimate M(6) G 0(61og(6) loglog(6)) from [3], 
this respectively gives the bounds 

0(d 2 n+1 Blog{B)) and 0(d 2 n+1 B \og 2 (B) loglog(B)). 

These estimates are far from optimal. A first important improvement (present in 
our code) consists in interpolating in the first place only the leading coefficients of 
the subresultants, and recover all other coefficients when needed. This is sufficient 
for the algorithms of Section [3J For instance, in the FFT case, the cost is reduced 
to 

0(d 2 n+1 B + d n+1 Blog(B)). 

Another desirable improvement would of course consist in using fast arithmetic based 
on Half-GCD techniques [IT] , with the goal of reducing the total cost to 0~(d n+ iB), 
which is the best known bound for computing the resultant, or a given subresultant. 
However, as of now, we do not have such a result, due to the possible splittings. 
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5.2 Solving two equations 



Our goal now is to estimate the cost of computing the polynomials R and G in 
the context of Formula [I] p. [6j We propose an approach where the computation of 
G essentially comes come free, once R has been computed. This is a substantial 
improvement compared to traditional methods, such as [TH [22] . which compute G 
without recycling the intermediate calculations of R. With the assumptions and 
notations of Section 15. 14 we saw that the resultant R can be computed in at most 
0(d n+ \B\og(B) + d^ +l B) operations in k. In many cases (typically, with random 
systems), G has degree one in y = x n+ \. Then, the GCD G can be computed within 
the same bound as the resultant. Besides, in this case, one can use the Half-GCD 
approach instead of computing all subresultants of P and Q. This leads to the 
following result in the bivariate case; we omit its proof here. 

Corollary 1 With n = 2, if V(H, P,Q) is empty and deg(G,y) = 1, then solving 
the input system P = Q = can be done in 0~(d\di) operations in k. 

5.3 Implementation of Regularize 

The operation Regularize specified in Section [2TT1 is a core routine in methods com- 
puting triangular decompositions. It has been used in the algorithms presented in 
Section HI Algorithms for this operation appear in [Tf|,[22]. 

The purpose of this section is to show how to realize efficiently this operation. For 
simplicity, we restrict ourselves to regular chains with zero-dimensional saturated 
ideals, in which case the separate operation of [H] and the regularize operation (22] 
are similar. For such a regular chain T in k[x] and a polynomial P G k[x] we 
denote by RegularizeDimO(P, T) the function call Regularize(P, T). In broad terms, 
it "separates" the points of V(T) that cancel P from those which do not. The output 
is a set of regular chains {T 1; . . . ,T e } such that the points of V(T) which cancel p 
are given by the Tj's modulo which p is null. 

Algorithm [1] differs from those with similar specification in [T^J [22] by the fact it 
creates opportunities for using modular methods and fast polynomial arithmetic. 
Our first trick is based on the following result (Theorem 1 in [1]): the polynomial p 
is invertible modulo T if and only if the iterated resultant of P with respect to T is 
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non-zero. The correctness of Algorithm [T] follows from this result, the specification of 
the operation RegularGcd and an inductive process. Similar proofs appear in [HJE2]- 
A proof and complexity analysis of Algorithm [1] will be reported in another article. 

The main novelty of Algorithm [1] is to employ the fast evaluation/interpolation 
strategy described in Section 15.11 In our implementation of Algorithm [fl at Step 
(6), we compute the "Scube" representing the subresultant chain of q and C v . This 
allows us to compute the resultant r and then to compute the regular GCDs (g, E) 
at Step (12) from the same "Scube". In this way, intermediate computations are 
recycled. Moreover, fast polynomial arithmetic is involved through the manipulation 
of the "Scube". 

Algorithm 1 

Input: T a normalized zero- dimensional regular chain and P a polynomial, both in 

k[Xi , • • • , %n\ ■ 

Output: See specification in Section POl 



RegularizeDimO(P,T) == 

(1) Results := 0; 

(2) for (q,C) G RegularizelnitDimO(P, T) do 



(3) if q G k then 

(4) Results := {C} U Results 

(5) else v :— mvar(g) 

(6) r ;= res(g, C v , v) 

(7) for D e RegularizeDimO(r, C <v ) do 

(8) s := Norma I Form (r, D) 

(9) if s ^ then 

(10) U := {DU {C v } U C >v } 

(11) Results := {U} U Results 

(12) else for (g,E) G RegularGcd (g, C v , D) do 

(13) g := NormalForm(#,£) 

(14) U :={EU{g}UD >v } 

(15) Results := {U} U Results 

(16) c := NormalForm(quo(C t; , g), E) 

(17) if deg(c,v) > then 
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(18) Results := RegularizeDimO(g, EUcU C>u) U Results 

(19) return Results 



In Algorithm [TJ a routine RegularizelnitialDimO is called, whose specification is given 
below. See [22] for an algorithm. 

Input: T a normalized zero- dimensional regular chain and p a polynomial, both in 
k[xi, . . . , x n \. 

Output: A set of pairs {{pi, Tj) | i — 1 • • • e}, in which pi is a polynomial and Tj is 
a regular chain, such that either pi is a constant or its initial is regular modulo 
sat(Tj), and p = Pi mod sat(Tj) holds; moreover we have T — > (T 1; . . . , T e ). 



6 Experimentation 

We have implemented in C language all the algorithms presented in the previous sec- 
tions. The corresponding functions rely on the asymptotically fast arithmetic opera- 
tions from our modpn library [TU|. For this new code, we have also realized a Maple 
interface, called FastArithmeticTools, which is a new module of the RegularChains 
library [18J. 

In this section, we compare the performance of our FastArithmeticTools com- 
mands with Maple's and Magma's existing counterparts. For Maple, we use its 
latest release, namely version 13; For Magma we use Version V2 . 15-4, which is the 
latest one at the time of writing this paper. However, for this release, the Magma 
commands TriangularDecomposition and Saturation appear to be some time 
much slower than in Version V2 . 14-8. When this happens, we provide timings for 
both versions. 

We have three test cases dealing respectively with the solving of bivariate systems, 
the solving of systems of two equations and the regularity testing of a polynomial 
w.r.t. a zerodimensional regular chain. In our experimentation all polynomial coef- 
ficients are in a prime field whose characteristic is a 30bit prime number. For each 
of our figure or table the "degree" is the total degree of any polynomial in the input 
system. All the benchmarks were conducted on a 64bit Intel Pentium VI Quad CPU 
2.40 GHZ machine with 4 MB cache and 3 GB main memory. 
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For the solving of bivariate systems we compare the command Triangularize of the 
RegularChains library to the command BivariateModularTriangularize of the 
module FastArithmeticTools. Indeed both commands have the same specification 
for such input systems. Note that Triangularize is a high-level generic code which 
applies to any type of input system and which does not rely on fast polynomial arith- 
metic or modular methods. On the contrary, BivariateModularTriangularize is 
specialized to bivariate systems (see Section 15.21 and Corollary [T]) is mainly im- 
plemented in C and is supported by the modpn library. BivariateModularTri- 
angularize is an instance of a more general fast algorithm called FastTrian- 
gularize; we use this second name in our figures. 

Since a triangular decomposition can be regarded as a "factored" lexicographic 
Grobner basis we also benchmark the computation of such bases in Maple and 
Magma. 
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Fi gUre 2: Generic dense bivariate systems. 



Figure [2] compares FastTriangularize and (lexicographic) Groebner : -Basis in 
Maple on generic dense input systems. On the largest input example the former 
solver is about 20 times faster than the latter. Figure [3] compares FastTrian- 
gularize and (lexicographic) Groebner : -Basis on highly non-equiprojectable dense 
input systems; for these systems the number of equiprojectable components is about 
half the degree of the variety. At the total degree 23 our solver is approximately 
100 times faster than Groebner : -Basis. Figure H] compares FastTriangularize, 
GroebnerBasis in MAGMA and TriangularDecomposition in MAGMA on the same 
set of highly non-equiprojectable dense input systems. Once again our solver out- 
performs its competitors. 
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Figure 4: Highly non-equiprojectable bivariate systems. 



For the solving of systems with two equations, we compare FastTriangularize (im- 
plementing in this case the algorithm described in Section I5T2]) with GroebnerBasis 
in Magma. On Figure [5] these two solvers are simply referred as Magma and 
Maple. For this benchmark the input systems are generic dense trivariate systems. 

Figures EJ [7J and M compare our fast regularity test algorithm (Algorithm [TJ) with the 
RegularChains library Regularize and its MAGMA counterpart. More precisely, in 
Magma, we first saturate the ideal generated by the input zerodimensional regular 
chain T with the input polynomial P using the Saturation command. Then the 
TriangularDecomposition command decomposes the output produced by the first 
step. The total degree of the input z-th polynomial in T is dj. For Figure [6] and 
Figure [7] the input T and P are random such that the intermediate computations do 
not split. In this "non-splitting" cases, our fast Regularize algorithm is significantly 
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Fi glire 5i Generic dense trivariate systems. 



faster than the other commands. 

For Figure [8]the input T and P are built such that many intermediate computations 
need to split. In this case, our fast Regularize algorithm is slightly slower than 
its MAGMA counterpart, but still much faster than the "generic" (non- modular and 
non-supported by modpn) Regularize command of the RegularChains library. The 
slow down w.r.t. the Magma code is due to the (large) overheads of the C - Maple 
interface, see [19] for details. 
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Figure 6: 2-variable random dense case. 
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Figure 7: 3-variablc random dense case. 
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Fig lire 8'. 3- variable dense case with many splittings. 



7 Conclusion 

The concept of a regular GCD extends the usual notion of polynomial GCD from 
polynomial rings over fields to polynomial rings modulo saturated ideals of regular 
chains. Regular GCDs play a central role in triangular decomposition methods. Tra- 
ditionally, regular GCDs are computed in a top-down manner, by adapting standard 
PRS techniques (Euclidean Algorithm, subresultant algorithms, . . . ). 

In this paper, we have examined the properties of regular GCDs of two polynomials 
w.r.t a regular chain. The theoretical results presented in Section [3] show that one 
can proceed in a bottom- up manner. This has three benefits described in Section [51 
First, this algorithm is well-suited to employ modular methods and fast polynomial 
arithmetic. Secondly, we avoid the repetition of (potentially expensive) intermediate 
computations. Lastly, we avoid, as much as possible, computing modulo regular 
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chains and use polynomial computations over the base field instead, while controlling 
expression swell. The experimental results reported in Section [6] illustrate the high 
efficiency of our algorithms. 
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